### Beginning ####
rm(list = ls())

# combine models from both studies that end up in the main text into a single table
pacman::p_load(readr, modelsummary)
mm_1 <- read_rds("./data_clean/study1_main_models_unscaled.rds")
mm_2 <- read_rds("./data_clean/study2_main_models_unscaled.rds")

psych_long_us <- read_dta("./data_clean/psych_clean_long.dta") %>% 
  rename(StoryTrue = true,
         # CRT = zCRT_all,
         CRT = CRT_all, #CRT IS NOW UNSCALED VERSION
         # AOT = AOT_all,
         AOT =  AOT_allUnscaled,
         # AOT_shortUnscaled,
         # AOT_longUnscaled,
         # RussianOrientation = russianOrientation,
         # ProRussiaPoliOrientation = zantiEuropeOrientation,
         RussianOrientation = XrussianOrientation,
         ProRussiaPoliOrientation = antiEuropeOrientation,
         Female = RSPSEX,
         Age = RSPAGE,
         Education = RSPEDUC,
         Income = income) %>% 
  mutate(Belief = response,
         AOT = (AOT - 1)/4)

psych_long2_us <- read_dta("./data_clean/psych_clean_long2.dta") %>% 
  mutate(RSPEDUC = factor(RSPEDUC),
         INCOM = factor(INCOM),
         oblast = as_factor(oblast),
         narr_topic = as_factor(narr_topic)) %>% 
  rename(StoryTrue = narr_true,
         # CRT = zCRT_all,
         CRT = crt,
         # AOT = AOT_all,
         AOT = AOT_allUnscaled,
         # RussianOrientation = russianOrientation,
         RussianOrientation = russianOrientation_unscaled,
         # RussiaNotAThreat = zthreat_russian,
         # RussiaNotAThreatUnscaled = RussiaNotAThreat,
         RussiaNotAThreat = RussiaThreat,
         Female = RSPSEX,
         Age = RSPAGE,
         Education = RSPEDUC,
         Income = INCOM,
         OblastStrata = strtcode) %>% 
  mutate(Belief = narr_believe,
         OblastStrata = as_factor(OblastStrata),
         CRT = CRT/4,
         AOT = (AOT - 1)/4)

vcov_mm_1 <- lapply(mm_1, cluster.vcov, ~ Respondent_Serial + question, parallel = TRUE)
vcov_mm_2 <- lapply(mm_2, cluster.vcov, ~ id + narrative, parallel = TRUE)

#doesn't allow you to combine variables, so try without.
cm <- c('(Intercept)' = 'Constant',
        'StoryTrue' = 'True',
        'CRT' = 'CRT',
        'AOT' = 'AOT',
        'RussianOrientation' = 'Russian Orientation',
        'ProRussiaPoliOrientation' = 'Russia Pol. Orientation',
        'CRT:StoryTrue' = 'CRT×True',
        'AOT:StoryTrue' = 'AOT×True',
        'CRT:RussianOrientation'= 'CRT×RussianOrientation',
        'AOT:RussianOrientation' = 'AOT×RussianOrientation',
        'StoryTrue:RussianOrientation' = 'True×RussianOrientation',
        'CRT:StoryTrue:RussianOrientation'    = 'CRT×True×RussianOrientation',
        'AOT:StoryTrue:RussianOrientation'    = 'AOT×True×RussianOrientation',
        'CRT:ProRussiaPoliOrientation' = 'CRT×Pro Russia Pol. Orientation',
        'AOT:ProRussiaPoliOrientation'  = 'AOT×Pro Russia Pol. Orientation',
       'StoryTrue:ProRussiaPoliOrientation' = 'True×Pro Russia Pol. Orientation',
       'CRT:StoryTrue:ProRussiaPoliOrientation' = 'CRT×True×Pro Russia Pol. Orientation',
       'AOT:StoryTrue:ProRussiaPoliOrientation' = 'AOT×True×Pro Russia Pol. Orientation')


cm2 <- c('(Intercept)' = 'Constant',
        'StoryTrue' = 'True',
        'CRT' = 'CRT',
        'AOT' = 'AOT',
        'RussianOrientation' = 'Russian Orientation',
        'ProRussiaPoliOrientation' = 'Russia Pol. Orientation',
        'RussiaNotAThreat' = 'Russia Not A Threat',
        'CRT:StoryTrue' = 'CRT×True',
        'AOT:StoryTrue' = 'AOT×True',
        'CRT:RussianOrientation'= 'CRT×RussianOrientation',
        'AOT:RussianOrientation' = 'AOT×RussianOrientation',
        'StoryTrue:RussianOrientation' = 'True×RussianOrientation',
        'CRT:StoryTrue:RussianOrientation'    = 'CRT×True×RussianOrientation',
        'AOT:StoryTrue:RussianOrientation'    = 'AOT×True×RussianOrientation',
        'CRT:RussiaNotAThreat' = 'CRT×Russia Not A Threat',
        'AOT:RussiaNotAThreat'  = 'AOT×Russia Not A Threat',
        'StoryTrue:RussiaNotAThreat' = 'True×Russia Not A Threat',
        'CRT:StoryTrue:RussiaNotAThreat' = 'CRT×True×Russia Not A Threat',
        'AOT:StoryTrue:RussiaNotAThreat' = 'AOT×True×Russia Not A Threat')

f <- function(x) format(round(x, 3), big.mark=",")
gm <- list(
  list("raw" = "nobs", "clean" = "Observations", "fmt" = f),
  list("raw" = "adj.r.squared", "clean" = "Adjusted R^2", "fmt" = f))

msummary(c(mm_1), output = "./tables/study_1_main_paper.html",
         title = "Study 1: Linear Models",
         notes	= c("Two-way standard errors clustered on Respondent and Narrative"),
         statistic_override = c(vcov_mm_1),
         coef_map = cm,
         gof_map = gm,
        escape = TRUE,
         stars = TRUE)   #tab_spanner(label = 'Literacy', columns = c('OLS 1', 'NBin 1')) %>% # add columns

names(mm_2) <- paste("Model", seq(1:length(mm_2)))

rows <- tibble::tribble(~term, ~`Model 1`, ~`Model 2`,~`Model 3`, ~`Model 4`,
          "Oblast Strata Fixed Effects", "Yes", "Yes", "Yes", "Yes")
attr(rows, 'position') <- c(37)
msummary(c(mm_2), 
         output = "./tables/study_2_main_paper.html",
         title = "Study 2: Linear Models",
         notes	= c("Two-way standard errors clustered on Respondent and Narrative"),
         statistic_override = c(vcov_mm_2),
         coef_map = cm2,
         gof_map = gm,
         escape = TRUE,
         add_rows = rows,
         stars = TRUE)   #tab_spanner(label = 'Literacy', columns = c('OLS 1', 'NBin 1')) %>% # add columns
